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Abstract 

Numerically optimized reduced descriptions of macromolecular liquids often present thermody- 
namic inconsistency with atomistic level descriptions even if the total correlation function, i.e. the 
structure, appears to be in agreement. An analytical expression for the effective potential between 
a pair of coarse-grained units is derived starting from the first-principles Ornstein-Zernike equation, 
for a polymer liquid where each chain is represented as a collection of interpenetrating blobs, with 
a variable number of blobs, n&, of size iV&. The potential is characterized by a long tail, slowly 
decaying with characteristic scaling exponent of This general result applies to any coarse- 

grained model of polymer melts with units larger than the persistence length, highlighting the 
importance of the long, repulsive, potential tail for the model to correctly predict both structural 
and thermodynamic properties of the macromolecular liquid. 
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Relevant structural and dynamical properties of polymer liquids take place over a large 
range of length (and time) scales. The computational requirements often render it imprac- 
tical or impossible to fully investigate such systems with atomistic level simulations. To 
overcome these limitations, many descriptions of polymeric liquids with reduced internal de- 
grees of freedom, or coarse-grained (CG) descriptions, have been proposed. The tremendous 
interest that CG methods have generated is due to their capability of speeding up simula 
tions, probing systems on larger length scales than conventional atomistic simulations, [l 

The basis for any coarse grained description is the specification of effective interaction 
potential energies or forces between coarse grained units. Despite the growing interest in 
CG methods and their proven successes, most CG models are still limited in their potential 
application because of the empirical character of the effective interaction potentials upon 
which they rely. While a few CG approaches have been formally derived, 9Ml2| most, for 
example the Iterative Boltzman Inversion (IBI) procedure where the mesoscale potential 
is optimized to reproduce the total correlation functions, and sometimes secondarily the 
pressure, jj] rely on numerical optimization of their CG parameters through comparison with 
the related atomistic simulations or with experimental data. Unfortunately, numerically 
optimized CG potentials are in principle neither transferable to different systems, nor to the 
same system in different thermodynamic conditions, and they can not assure consistency for 
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properties different from the ones against which their parameters have been optimized 
This limits dramatically their generality and convenience. 

In this letter we present an analytical formalism for the potential describing the interac- 
tion between CG units, for a model of a polymer liquid where each molecule is represented 
as a chain of interpenetrating spheres with variable size and number, in variable thermody- 
namic conditions: we show that the CG model presented and the related analytical potential, 
ensure structural and thermodynamic consistency with the atomistic description. The model 
is general, as the potential is fully transferable and it is applicable to liquids of polymers with 
different molecular structure. The analytical form of the potential allows for the character- 
ization of some general features of how thermodynamic properties and structure depend on 
the shape of the intermolecular potential between CG units. 

In a polymeric liquid of monomer density, p m , and number of polymers n, the size of 
the polymeric chain is defined by its radius-of-gyration, R 2 g = Nl 2 /6, with I the effective 
segment length between the center-of-mass (corns) of two monomers, and N the number 
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of monomers in a chain. Here we are concerned about liquids of polyethylene for which 
I = 0.152 nm. Each polymer is represented as a chain of superimposing soft blobs, with 
Nb the number of monomers in one blob, rib = iV/iVj, the number of blobs in one chain 
and blob size characterized by, R g b = R g j ^/nb, The correlation functions that describe the 
static structure of the blob CG model have been derived from the formal solution of a 
generalized Ornstein-Zernike equation where monomeric sites are assumed to be real sites, 
and CG sites are assumed to be fictitious sites, in an extension of Krakoviack's et al. orig- 



inal approach, Il5j and including the thread model polymer reference interaction site 

model (PRISM) monomeric description. jl^l . \y\ By assuming a Gaussian distribution of the 

monomers inside a blob and monomer interactions much shorter ranged than the size of the 

blob, which are reasonable approximation for subchains with Nb larger than the persistence 

length (for polyethylene JV& > 30), an analytical blob-blob total correlation function, h bb (r), 

was derived, and shown to be in quantitative agreement with the total correlation function 

□ 

from atomistic simulations. |17l] 

Starting from h bb {r) and given that the monomer direct correlation function, c mm (r), has 
a range shorter than the size of the CG unit, an analytical solution of the CG potential 
is derived by approximating the Fourier transform of c mm (r) by its zero wavevector value, 
Co = 47r J °° r 2 c mm (V). The effective direct correlation between blobs is given in Fourier space 
by 

c bb (k) = NbTb fe(*0/"£(*0) 2 m 

where Tb = —Nbp m co where cq < 0, and Cl b ™(k) is the intramolecular correlation function be- 
tween the corns of a blob and a monomer, averaged over the chain and normalized such that 
its value at k = is equal to one. The sharp peak around k = leads to the enhanced sensi- 
tivity of the function to the quality of the approximations used in the chain model, and made 
numerical solutions of the potential necessary in our previous work. The total correlation 
function in this blob description gives isothermal compressibility, kj- = [ksT 'p m (l+I \) / 'N]^ 1 , 
with ks the Boltzmann constant and T the temperature, which is consistent with the com- 



pressibility in the atomistic description[16[j. Intra- and inter-molecular pair distribution 
functions and the structure factor, are consistent with their atomistic counterpart for dis- 
tances r > R g b and wave vectors k < 2ir /R g b- 

Assuming that a coarse grained unit contains a number of monomers sufficiently large to 
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follow a Gaussian space distribution, and that the density and the interaction strength are 
large enough that the product Tf, » 1, the contribution to the inverse transform integral for 
large wavevectors (k » 1/R g ^) is negligible for r > R g t,. In this limit, the direct correlation 
for k << 1/Rgb has the simple rational function limiting form c bb (k) oc 1/(1 + TR^k 4 ) as 
r& — y oo. Approximating the effective direct correlation by this form for all wavevectors 
introduces very little error, allowing for a simple approximation for the functional form in 
real space. The accuracy at intermediate T values can be improved by taking this as the 
zeroth order term of an asymptotic expansion in 1/y/Tl about Tb — > oo. 

The effective potential is then derived by applying the approximation V bb (r) ps 
-k B T[c bb (r) - h bb (r) + ln(l + h bb (r))] » -k B Tc bb (r), the limiting form of the Hyper-Netted 
Chain approximation valid when \h bb (r)\ << 1 everywhere. This approximation holds for 



For r > R 
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soft potentials in the limit of high densities and long chains of interest here, 
the intermolecular blob potential, which is the needed input for the mesoscale simulations 
of the coarse-grained polymer liquid, is given by 

V bb (r) ~ kr>T\ ( 45 ^ r ^ 4 ^ sinjQr) Qr 

V [r) ~ K B 1 [ y^^tfj Qr e 
VMb i [/iqn3 



672 7 rp m r? 7T J R3 



[{UQUQr - 4))cos(Qr) (2) 



b " 9 6 / 

Qr 



945 + 1 / 4 9r3 rsm{Qr) + ^_cos{Qr)\- 



r l/4 I ^w. v ^, , , i /4 ^^. n Qr 

b / b 



where Q = 5 1 / 4 y / 3/2/ (I^) 1 / 4 , Q r s = 5 1 / 4 y / 3/2, and r is in units of the blob radius-of- 
gyration, _R 9 &. When rib = 1 Eq-© represents the effective potential between the center-of- 
mass of two polymers in a melt. 

Because the potential is formally expressed as a function of the molecular and thermody- 
namic parameters, Eq.© is in principle general and applicable to polymer melts in different 
thermodynamic conditions and with diverse macromolecular structures. 

The range of the effective potentials between CG units scales beyond the effective blob 
radius of gyration, R g b, and decays with the number of monomers per blob as N b 1/4 . The 
observed scaling behavior can be explained by considering that in a liquid the total effective 
correlation between two sites (atomistic or CG), and its related potential, can be regarded as 
"propagating" through sequences of direct pair interactions following the Ornstein-Zernike 
integral equation theory. These many-body contributions to the pair interaction are not 
simply additive, and once mapped into the OZ pair interaction, they result in a slowly de- 
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caying tail. Because of the Gaussian statistics that applies to the structure of long polymeric 
chains, the interaction between any intermolecular pair of blobs statistically propagates in 
the shared volume as a random walk following the random path of effective CG sites. Given 
that in the relevant volume Vj oc R^ b oc N b ^ 2 l 3 , the number of effective CG sites is of the 
order of n' b oc p m R z gb /N b oc N^pml 3 , which leads to the observed scaling of iV b 1//4 . 

From the analytical form of the effective potential, Eq.(j2]), the pressure of the system 
can be calculated explicitly. Specifically in the high density, long chain limit, the pressure 
calculated in the virial route reduces to the simple expression 

P , _ Nc p m 

Pc k B T~ l 2 • {6} 
_Qxn/N is the chain density. Eq.([3]) is in agreement with the monomer level 



where p c 



description, 16j and does not depend on the level of coarse-graining of the model selected. It 



should be noted that because both the pressure and the compressibility are consistent with 
monomer level PRISM, the equations of state predicted from the virial and compressibility 
routes will show the same small inconsistencies between the routes, which are inherent to the 
approximations in integral equation theory. 
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The important point is that the pressure 



and compressibility are left unchanged from the monomer level description by the coarse 
graining procedure on all levels. 

To test the self consistency and the predictions of the effective potentials we performed 
molecular dynamics (MD) simulations of polymer melts using the LAMMPS simulation 
package, 19| in parallel on the SDSC Trestles cluster accessed through the XSEDE project. 
Simulations were performed both at the atomistic level, by adopting the traditional united 
atom model (UA-MD) with the established set of UA potentials (intermolecular Lennard- 
26| harmonic bonds and angle potentials) and at the coarse-grained mesoscale 
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level (MS-MD) where each chain was represented as a collection of soft colloidal particles. 
MS-MD simulations for the n b = 1 case have been described previously 2CH22l| . To extend 
the model to cases where n b > 1, the effective bond between adjacent blobs was taken to be 
Vbond = SksTr 2 / (8R^ b ) + V bb (r) + ksTln^l + h bb (r)) , which enforces the correct distributions 



between adjacent blobs[l3], and an angle potential between sequential triples on each chain 
which likewise enforces the correct angular probability distribution. 23| Blobs more than two 
apart on each chain were taken to interact via the intermolecular pair potential V bb (r) of 
Eq.([2]). All Mesoscale simulations used here are performed in the NVE (microcanonical) 



ensemble. Values of the Co parameter, entering 1^ in Eq.([2]) and the MS-MD, were taken 
from the UA simulations for short chains. For systems with large N, which were too slow to 
relax to be accessible through UA-MD, the Cq parameter was extrapolated from the available 
UA data at small N at the same monomer density, using the form predicted from monomer 
level numerical PRISM theory, Co = a + b/N, with a and b optimized parameters. The values 
of Co obtained from this procedure were found to be generally consistent with calculations 
using long-established numerical monomer level PRISM models when an attractive part to 
the monomer potential is included, [lij 

Figure [T] displays the effective intermolecular force (— dV bb /dr) between blobs, comparing 
the numerical solution of the potential with the approximate expression, Eq.(J2]). The ap- 
proximate expression represents correctly the force in real space going from the peak to the 
long tail, which is the essential information needed for the correct calculation of thermody- 
namic properties. The inset in Figure [1] shows the normalized Virial density (N~ 2 r 3 F bb (r)), 
whose integral is proportional to the pressure, and which can be seen to be dominated by 
the tail region of the potential (r > R g b). The peak of the force decreases with increasing 
blob length, 7V&, at constant density, p m , but its range increases in such a way that the Virial 
integral ultimately reaches a plateau, which corresponds to the leveling off of the pressure, 
as shown in Figure |2j The inset of Figure [1] also shows that the effective potential has a 
small attractive part at long range: such an attractive component is a necessary condition 
for a stable liquid to form, as a gas phase would be the state of lowest free energy at any 
temperature for a system of purely repulsive particles with no additional constraints. This 
attractive contribution is partially of entropic origin and partially due to the attractive 
component of the inter- monomer potential. Figure [2] shows that all simulations of the same 
system, independent of the level of coarse graining, generate not only consistent structure, 
as seen before, but also consistent pressure. 

The scaling of the potential is a property of the representation; increasing the number 
of monomers in each coarse grained unit on each chain, with a fixed system density and 
monomer interaction, results in an effective potential that grows long-ranged because it 
captures the average effect of the many-polymer correlations mapped into the effective pair 
interactions. This is a consequence of the fractal dimension of polymers, which is of order 
two, as blob volume increases with Nb faster than the chain can fill it, leading to an increasing 
with Nb of the level of chain interpenetration. If however, melts of increasingly long chains 
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FIG. 1: Force curves for the rib = 1 model (solid lines) and the approximated Eq.([2]) (squares) 
for (top to bottom) N = 100 (red), N = 200 (blue), N = 500 (green), and N = 1000 (brown), 
along with force curves for models with a fixed blob length, 7V& = 50, and increasing numbers of 
blob (dashed lines): rib = 2 (red), = 4 (blue), rib = 10 (green), and rib = 20 (brown). Inset: 
Normalized Virial density (N~ 2 r s F bb (r)) for the numerical forces and their analytical forms for 
rib = 1- All systems are at 400K and density p m = 0.03355A" 3 . 

at the same density are represented by increasing numbers of soft blobs of fixed iV&, the 
range of the potential between blobs limits to a fixed value with increasing chain length as 
the direct correlation parameter Co approaches its limit. 

Analysis of these results shows that the range of the effective potential is highly density 
dependent, both explicitly and indirectly through the direct correlation parameter, cq, as 
shown in Figure |3j Also reported for comparison is the force generated from the potential 
of mean force, which is a very poor representation of the "real" force at high densities. 

Finally Figure H] compares predictions of pressure as a function of density for two samples 
with increasing chain length, i.e. iV = 44 and N = 100. Data from MS-MD of the CG 
soft-blob representation show excellent agreement with data from UA-MD in the range 
where UA-MD were performed. The agreement appears to be independent of the level of 
CG representation that is adopted, as soft-spheres and chains of soft-blobs have consistent 
pressure across the different levels of coarse graining. The calculations are dominated by 
the presence of the long-ranged tail of the potential, further validating the proposed CG 
description. The range of separations over which the effective potential must represent 
the average of many-polymer effects increases dramatically with density. While all levels 
of representation accurately reflect whole-system thermodynamic averages and structural 
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FIG. 2: Pressure measured from simulations with input CG forces from Figure Q] where = 1 
(orange squares connected by dashed guide lines), and Nf, = 50 with rib > 1 (open blue circles). 
Also depicted are UA-MD simulations (green X symbols) for systems which relax fast enough for 
UA-MD to be feasible (N < 200). Data are collected for increasing degree of polymerization 
N, at the same values of N as Fig.(pQ). All simulations were performed at 400K and density 

p m = 0.03355A- 3 . 
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FIG. 3: Numerical force curves for iV = 100, n& = 1 models of PE with increasing densities 
(bottom to top) 0.032261" 3 (red), 0.033551" 3 (blue), 0.03441 A" 3 (green), 0.03656A" 3 (black), 
and 0. 03871 A~ 3 (gray), along with corresponding mean force curves (dashed lines). Insert: detail 
of the attractive contribution to the force (densities increase left to right). All systems are at 400K 



pair correlations, a coarse grained description may be unable to resolve processes below the 
length scale of the potential tail, due to this averaging effect. 

The difference between the scaling predictions for the effective potential and the potential 
of mean force also has very important implications for theories of polymer melts. The 
strength of the potential of mean force, w(r) = — /csT/n(l + h bb (r)), is found to scale for 
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FIG. 4: Pressures measured in MS-MD simulations versus density, for N = 44 (triangles) and N = 
100 (circles), performed with effective potentials where cq is determined from UA-MD simulations. 
Also reported are values from the corresponding UA-MD simulations (crosses for N = 44 and stars 
for N = 100). All simulations were performed at 400K Connecting lines are added as a guide to 
the eye. 

long chains, or large blobs, and high densities as 1/ '(p m yN) and its range as N 1 ' 2 with 
Co becoming irrelevant and w(r) vanishing in the infinite chain limit. Use of the w(r) in 
thermodynamic relations therefore results in vanishing energies and pressures in the long 
chain/high density limit, which would imply the irrelevance of intermolecular interactions 
and allow reduction of the system to a single chain problem. While the effective pair potential 
at contact, V bb (0), still vanishes in the infinite chain limit, its tail increases in such a way that 
the Virial integral does not vanish. Likewise the average effective energy in the system also 
does not vanish, and furthermore both depend on the monomer level direct correlation, and 
therefore on the monomer interaction potential. This implies that even in the infinite chain 



limit, intermolecular potentials do not become irrelevant, and therefore cannot general 



neglected, as it is conventionally done in the description of polymer melt dynamics. 23] 



y be 



In this letter we have presented an analytical characterization of the soft pair potentials 
that arise in high level coarse grained descriptions of interpenetrating polymer melts, i.e. 
for soft-sphere and chains of soft-blobs coarse graining mappings of polymer melts. The 
form of the potential and the pressure are obtained directly in terms of parameters from the 
monomer level theory, without a need for phenomenological correction forms in the poten- 
tial and numerical optimization procedures. The effective potentials for these CG models 
show characteristic long-ranged "tails" that scale non-trivially with chain length, density 
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and monomer interaction strength. This family of CG models guarantees the consistency 
of the structure and thermodynamics of the macromolecular liquid in any level of soft rep- 
resentation, which is relevant in the modeling of complex polymeric liquids, as well as in 
the design of reliable multiscale modeling approaches to capture relevant phenomena that 
may occur on many different length scales. Dynamical properties will be accelerated in 
the coarse-grained representation, and will have to be properly rescaled to reconstruct the 



realistic dynamics. 



22] 



We acknowledge support from the National Science Foundation through grant DMR- 
0804145. Computational resources were provided by Trestles through the XSEDE project 
supported by NSF. This research was supported in part by the National Science Foundation 
under Grant No. PHY1 1-25915. 



[1] L. Jacobson and V. Molinero, J. Phys. Chem. B 114, 7302 (2010). 
[2] J. J. de Pablo Ann. Rev. Phys. Chem. 62, 555 (2011). 

[3] A. Savelyev, and G. A. Papoian, Proc. Natl. Acad. Sci. USA 107, 20345 (2010). 
[4] P. Carbone, H. A. Karimi Varzaneh, X.Y. Chen and F. Miiller-Plathe, J. Chem. Phys. 128, 
0649041 (2008). 

[5] V. A. Harmandaris, N. P. Adhikari, N. F. A. van der Vegt, and K. Kremer Macromol. 39, 
6708 (2006). 

[6] B. G. Levine, D. N. Le Bard, R. De Vane, W. Shinoda, A. Kohlmeyer, and M. L. Kein, J. 

Chem. Theory Comp. 7, 4135 (2011). 
[7] T. Vettorel, G. Besold, and K. Kremer, Soft Matter 6, 2282 (2009). 
[8] G. DAdamo, A. Pelissetto, and C. Pierleoni, Soft Matter 8, 5151 (2012). 
[9] S. Izvekov, and G. A. Voth J. Phys. Chem. B 109, 2469 (2005). 
[10] J. F. Rudzinski and W. G. Noid J. Chem. Phys. 135 214101 (2011). 

[11] P. G. Bolhuis, A. A. Louis, J.-P. Hansen, and E. J. Meijer. J. Chem. Phys.114, 4296 (2001). 
[12] G. Yatsenko, E. J. Sambriski, M. A. Nemirovskaya, and M. Guenza, Phys. Rev. Lett. 93, 
257803 (2004). 

[13] A. A. Louis J. Phys.: Condens. Matter 14, 9187 (2002). 

[14] V. Krakoviack, J.-P. Hansen, and A. A. Louis, Europhys. Lett. 58, 53 (2002). 



10 



[15] V. Krakoviack, B. Rotenberg, J. -P. Hansen J. Phys. Chem. B 108, 6697 (2004). 

[16] K. S. Schweizer and J. G. Curro, Chemical Physics 149, 105 (1990). 

[17] A. J. Clark, and M. G. Guenza J. Chem. Phys. 132, 044902 (2010). 

[18] J. -P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic, New York) 2 Ed. 
(1990). 

[19] S. Plimpton J. Comput. Phys. 117, 1 (1995). 

[20] J. McCarty, I. Y. Lyubimov, and M. G. Guenza J. Phys. Chem. B 113, 11876 (2009). 

[21] J. McCarty, I. Y. Lyubimov, and M. G. Guenza Macromol. 43, 3964 (2010). 

[22] I. Y. Lyubimov, and M. G. Guenza Phys. Rev. E 84, 031801 (2011). 

[23] M. Laso, H. C. Ottinger and U. W. Suter J. Chem. Phys. 95 2178 (1991). 

[24] E. K. Watkins and W. L. Jorgensen, J. Phys. Chem. A 105 4118 (2001). 

[25] M. Putz, J. G. Curro, and G. S. Grest, J. Chem. Phys. 114 2847 (2001). 

[26] M. G. Martin and J. I. Siepmann, J. Phys. Chem. B 102 2569 (1998). 

[27] M. Doi, and S. F. Edwards, The Theory of Polymer Dynamics (Oxford University, Oxford, 
1986). 



11 



